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ABSTRACT 


Previously we showed that a substantially misaligned viscous accretion disk 
with pressure that orbits around one component of a binary system can undergo 
global damped Kozai-Lidov (KL) oscillations. These oscillations produce peri¬ 
odic exchanges of the disk eccentricity with inclination. The disk KL mechanism 
is quite robust and operates over a wide range of binary and disk parameters. 
However, the effects of self-gravity, which are expected to suppress the KL oscil¬ 
lations for sufficiently massive disks, were ignored. Here, we analyze the effects of 
disk self-gravity by means of hydrodynamic simulations and compare the results 
with the expectations of analytic theory. The disk mass required for suppres¬ 
sion in the simulations is a few percent of the mass of the central star and this 
roughly agrees with an analytical estimate. The conditions for suppression of the 
KL oscillations in the simulations are close to requiring that the disk be gravita¬ 
tionally unstable. We discuss some implications of our results for the dynamics 
of protoplanetary disks and the related planet formation. 

Subject headings: accretion, accretion disks - binaries; general - hydrodynamics - 
planets and satellites: formation 
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Introduction 


Estimates sug gest that 40% to 50% of observed exoplanets are in binary systems 


flHorch et al 


2014i) . As the birthplace of the planets, protoplanetary disks are fnndamental 


to explaining planet formation and thus it is important to understand their evolution 
in a binary system. Recent ALMA observations revealed large mutual inclinations 


(greater than 60°' 


Jensen fc Akeson 


between the circumstellar disks around binary system components (e.g., 


2014J: 


Williams et ah 


2014i) . Although the binary orbital planes in these 


systems have not yet been identihed, it is possible that at least one of the disks in each 
system is signihcantly inclined ( > 45°) with respect to the binary orbit. 

The Kozai-Lidov (KL) mechanism is a well known process that occurs when a ballistic 
object orbits around one component of a binary system and the object’s orbital plane is 


antially misa ligned {i > = 39°) with respect to the binary orbital plane flKozai 


subst 

antiallv misaligned (i 

1962; 

Lidov 

1962 

). In this 


19621) . In this process, the orbital eccentricity and inclination of the object 


undergo periodical exchange. An object on an initially circular orbit periodically gains and 


loses eccentricity. There have be en a large number of works on this men 


hrst discovery in the 1960s (e.g., 


1998 : 


Ford et al. 


2000 


Holman et a. 


Lithwick fc Naoz 


2011 


1997 


Innanen et al. 


Naoz et al. 


2013al: 


1997 


lanism since its 


Liu et al. 


K^iseleva et al. 


2015 


Antognini 


most 

notablv here on the formation of hot 

upiter planets (e.g.. Wu & Murray 

200 

3; 

Takeda & Rasio 2005; Fabrvckv & Tremaine 

2007: Naoz et al. 2013b; Petrovich 

2015; 

Rice 

2015 

)• 


While the KL cycle of orbiting objects has been extensively studied, only very recently 


was it found t o operate on a fluic 


simulations, (IMartin et al. 


2014 


disk with pressure and viscosity by means of hydrodynamic 
, hereafter Paper I). Thus, an initially circular disk can 


become highly eccentric if its initial tilt exceeds about 45° with respect to the binary orbital 
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plane flFu et al.ll2015l. hereafter Paper II). Due to the efficient radial communication via 
either disk pressure or viscosity, the disk KL cycle for a typical protostellar disk operates 
in a global, coherent fashion, and the disk remains quite flat (unwarped) throughout the 
process. The eccentricity growth is fairly uniform in radius across the disk. Unlike the 
ballistic object case, the disk oscillations are damped due to viscous dissipation, likely 
involving shocks in the disk. When the oscillations fully damp, the disk is tilted at the 
critical KL angle (Ur ~ 40°) or somewhat less and the disk is approximately circular. After 


this stage, the circular disk evolves towards alignment with the bi nary orbita . 


tidal dissipation effects associated with turbulent viscosity (e.g.. 


King et ah 


plan e through 


2 OI 3 I) . Paper 


II extended the work of Paper I by surveying a large range of binary and disk parameter 
space, including the initial disk inclination, temperature, viscosity, size, surface density 
profile, and the binary mass ratio and eccentricity. The disk KL cycle is a fairly robust 
process that can occur under many different binary-disk conditions. However, Paper II did 
not consider the effects of self-gravity. 

At early times in the protoplanetary disk evolution, the disk mass is large, several 
percent of the mass of the central star, and the disk may experience dynamical effects of 
self-gravity. Depending on the properties of the disk, it is possible that self-gravity may be 


sufficiently strong to prevent KL oscillations from operating. Instea d, the disk wif 


Batvgin 

2012; 

Lai 

0 

_ 


ar obj ect (as it does for lower inclination disks) flBatvgin et al 


globa lly 


2011 


20141) . This suppression is expected when the disk apsidal precession 


rate due to self-gravity dominates over the apsidal precession rate due to the gravitational 
effects of the binary that cause the KL oscillations. The suppression may not continue as 
the disk evolves. As material drains on to the central star, the self-gravity weakens and the 
KL cycles can begin if the disk remains sufficiently misaligned. 


Differential precession due to the binary acts to disrupt a disk. For a disk to process 
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coherently, internal torqnes are reqni red to operate 


the di fferential precession timescale flLarwood et ah 


globa 


ly and commnnicate faster tha n 


19961) . The disk model of 


Batygin 


(120121) maintained disk coherence by means of disk self-gr avity and omitted the effects of 


pressure and viscosity. 


Xiang-Gruess fc Papaloizoul (120141 1 pointed out that pressure alone 


allows a disk to process coherently and thus self-gravity is not necessarily required. In 
this paper, for the hrst time we present results on the long term evolution of misaligned 
self-gravitating disks that have pressure and viscosity. 


2. Numerical Simulations 

In this section, we describe three-dimensional simulations of an inclined fluid disk 


We use the smoot 

K- 

led particle hydrodynamic fSPH) code PHANTOM ( 

Lodato & Price 

2010 ; 

Price & Federrath 

2010 ; 

Nixon 

2012 ; 

Nixon & King 

2012 ; 

Price 

2012 

Nixon et ah 

2013) 


The simulation setup is similar to that in Paper 11. An initially circular disk orbits around 
a central binary component of mass Me under the influence of the perturber binary 
component of mass Mp. The disk and stars interact gravitationally and respond to these 
interactions. In particular, the orbit of the binary is affected by its gravitational interactions 
with the disk, although the changes to the binary orbit are small. The orbital plane of the 
disk is initially inclined to the binary orbital plane by 50°. We adopt a locally isothermal 
equation of state and an explicit accretion disk viscosity. We include a nonlinear term with 
a coefficient /3 av = 2 (AV stands for artiheial viscosity) in order to suppress interparticle 
penetration, as is standard in SPH codes. The disk sound spe ed is oc and the in itial 


surface density distribution is S oc r such that both a (IShakura fc Sunvaev 


and t he smoothing length (h) /H are constant over the disk radius, r (iLodato &: Pringle 


1973) 


20071) . We start the simulations with 1 x 10® SPH particles in the disk. The PHANTOM code 

























































adopts a cubic spline kernel as the smoothing kernel. The number of neighbors is roughly 
constant at A^neigh ~ 58. The disk initially extends from radius rjn to radius rout- The inner 
boundary of the simulated region is set to the initial disk inner disk radius rju. As particles 
move to r < rjn, they are removed from the simulation while their mass and momentum are 
added to the central star. We also impose an inner boundary radius around the perturbing 
companion, since some disk mass can be transferred to that component. In terms of the 
binary orbital size Ob, the initial disk inner and outer radii in our simulations are located at 
Tin = 0.025ab and rout = 0.25ab, respectively. The initial outer disk radius is chosen to be 


the tidal truncation radius of a coplanar disk llPaczvnski 


19771). Howeve r, because the tidal 


2015fl . the disk initially 


torques on a misaligned disk are weaker flLubow. Martin fc Nixon 
expands somewhat beyond this radius. 

We describe simulations with two sets of parameters for the disk aspect ratio and 
viscosity. The hrst is H/R = 0.1 and a = 0.01 and the second is H/R = 0.06 and 
a = 0.01 {H/R is evaluated at disk inner edge rin). The disk is resolved with shell-averaged 
smoothing length per scale height (h) /H k. 0.26 and 0.37, respectively. For each set of 
H/R and a, we consider three different initial disk masses Ma = O.OOIM, O.OIM and 
0.03M (in units of the total binary mass M = Me + Afp). In Papers I and II, we ignored 
disk self-gravity given the low disk mass used there (O.OOIM). In this paper, we take into 
account the effect of disk self-gravity and compare simulations with and without self-gravity 


for the same disk mass. The algorithm for 


PHANTOM is described in 


he S PH implementation of self-gravity in 


Price fc Monaghan! (120071) . which discusses numerical tests based 


on the radial oscillations and the static structure of a polytrope. The self-gravity algorithm 


of PHANTOM has also been used to stuc 


20081 ) and star clusters (jPrice fc Bate 


y the formation of giant molecular clouds flDobbs 


20081) . 


Figure [T] shows the evolution of the disk eccentricity (upper row), inclination (middle 





















Table 1. Parameters of the SPH simulations for equal mass binary systems with total 

mass of M and separation of Ob 


Binary and Disk parameters Symbol Value 


Mass of each binary component 
Binary orbital eccentricity 
Initial number of particles 
Initial disk mass 
Initial disk outer radius 
Initial disk inner radius 
Mass accretion radius 
Disk viscosity parameter 
Disk aspect ratio 

Initial disk surface density profile S oc r 
Initial disk inclination 


Mp/M = Mc/M 

0.5 

Cb 

0 

N 

10® 

-Mdisk/ M 

[0.001, 0.01, 0.03] 

rout/^b 

0.25 

rin/ah 

0.025 

race/ O-b 

0.025 

a 

[0.01, 0.1] 

H/r (r = Tin) 

[0.06, 0.1] 

7 

1.5 

*0 

50° 






No Self-gravity With Self-gravity 



Fig. 1. — Evolution of the eccentricity (upper row), inclination (middle row), and total mass 
(bottom row) of the circumprimary disk. The left column is for simulations without disk self¬ 
gravity, whereas the right column is for simulations that take into account disk self-gravity. 
The eccentricity and inclination are measured at a radius of r = 0.2ab (the initial disk radial 
range is 0.025ab <r< 0.25ab), where Ob is the binary semi-major axis and Pb is the binary 
orbital period. Within each plot, each line represents a simulation with a different initial 
disk mass. The units of mass are that of the total binary mass, M. All lines are averaged 
over one binary orbital period. The initial number of SPH particles is 10®. The disk aspect 
ratio at the inner edge is H/R = 0.1. The disk viscosity parameter is a = 0.01 and the 
initial disk inclination is i = 50°. 
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row), and mass (bottom row) for three different initial disk masses. The eccentricity and 
inclination are measnred at a representative radins r = 0.2ab, since the disk remains flat 
and the eccentricity is ronghly constant with radins. The left and right colnmns show 
simnlations withont and with disk self-gravity, respectively. In all cases, disk gravity acts 
on the stars and slightly affects the binary orbit. In post-processing the simnlation ontpnts, 
we compnte the orbital eccentricity and inclination for each particle nsing its position and 
velocity information, then divide the disk into 100 radial bins and calcnlate the mean 
properties of the particles within each bin. Resnlts for the eccentricity and inclination are 
taken from the radial bin centered on radins 0.2ab. This is similar to the method nsed in 
Papers I and II, except now we average the properties over one binary orbit to smooth 
ont the flnctnations on the binary orbital timescale. Note that the initial disk eccentricity 
is non-zero (e ~ 0.07) in the npper row, even thongh the disk is initially set np to be 
circnlar. This non-zero eccentricity valne is a conseqnence of the way we calcnlate the 
particle eccentricity and inclination (Eqnations [7] and [8] in Paper II) that treats particles 
as ballistic, i.e., assnming particles only feel the gravitational force of the central object, 
whereas particles actnally are also inflnenced by the disk pressnre force. Therefore, the 
non-zero initial disk eccentricities we see here are pnrely an artifact of onr ignoring the 
disk pressnre in the orbital elements calcnlation. The effect is more prominent here than 
in the simnlations shown in Papers I and II becanse here the disk has higher temperatnre, 
H/R = 0.1 (compared with H/R = 0.035). However, it is still small enongh not to affect 
onr general interpretation, since we are focnsed on the major changes of the disk properties 
dnring KL oscillations. 

As seen in the left colnmn of FignrelH where disk self-gravity is not inclnded, increasing 
the disk mass does not signihcantly affect the disk orbital elements dnring the KL cycle. 
The period and amplitnde of the oscillations are similar, althongh there is some variation 
in the time of the hrst peak. The difference in the onset of the initial KL cycle is dne to 
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the fact that disk mass slightly affects the companion’s orbit. The timing of the onset is 
sensitive to the initial conditions. Even with = 0.03 M, the peak eccentricity valne of 
e ~ 0.42 is almost the same as that for the simnlations with = 0.001 M and 0.01 M 
with a peak e = 0.44. As seen in the right colnmn of the fignre, for the lower mass disks, 
self-gravity generally has little effect on disk eccentricity and inclination. However, the 
effects of the KL cycle are significantly rednced at the highest disk mass of = 0.03M , 
which has an initial disk minimnm Toomre Q ~ 2.6. The peak eccentricity in that case is 
e ~ 0.2 compared to e ~ 0.42 with Md = O.OOIM and O.OIM. The bottom row shows that 
by the end of the simnlation at time t = 29P]^, the disk in general has lost more than half of 
its initial mass. The majority of the lost mass is accreted on to the central object. A small 
fraction of the mass is ejected from the disk, mostly ending np aronnd the companion binary 
star. Comparing simnlations with the same initial disk mass, the rnn with self-gravity has 
more mass remaining at the end of the simulation than the one without self-gravity. This 
effect may be a consequence of the disk self-gravity helping to retain particles, and also 
reducing the particle loss at inner disk boundary, since the KL cycle is weakened. 

Figure |2] shows edge-on views of the disk in the two simulations with Md = 0.03 M, one 
with and the other without disk self-gravity. The left panel shows the initial conditions that 
were applied to both simulations. The middle and right panels are edge-on views after the 
disks have undergone nodal precession by 180°. The middle panel is from the run without 
disk self-gravity (the blue curves in the left column of Figure [T]) at a time of t = 14.9Pb, 
while the right panel is from the run that includes disk self-gravity (the blue curves in the 
right column of Figure [T]) at a time of t = 15.4Pb- These runs have different disk nodal 
precession rates and as a result the middle and right panels are shown at slightly different 
epochs. Without self-gravity, the middle panel shows a substantially lopsided disk due to 
the KL-driven eccentricity growth. In the right panel, with disk self-gravity, the disk is still 
eccentric, but not as eccentric as in the middle panel. Although at an earlier time {t ~ 7Pb) 
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x-z plane ^_q p 

\ 

I_I without SG with SG 

0.2 at _; 

Fig. 2.— View of the disk towards the x — z plane for the two simulations with Md = 0.03 M, 
H/R = 0.1 and a = 0.01, with and without self-gravity. The binary orbit is in the x — y 
plane (i.e., the perturbing object moves into and out of the page). The left panel shows 
the initial disk setup that is the same for both simulations. The middle (right) panel shows 
the disk without (with) self-gravity after it has undergone nodal precession by 180°. The 
central mass is denoted by the black dot. The color coding is for the logarithm base 10 of 
the column density (i.e., density integrated along the line of sight) in units of {Mc + M^)/a^. 
(Color online) 
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the disk has a slightly higher eccentricity, the average disk eccentricity is clearly lower when 
self-gravity is included (see comparison of the blue curves in the top two panels of Figure 
H]). The main point here is that in both panels, the disk remains fairly flat (no signihcant 
warping) due to the efficient radial communication (provided mainly by pressure in these 
cases; see discussions in Paper II). 

Figure [3] is similar to Figured] except for a different disk aspect ratio, H/R = 0.06, and 
viscosity, a = 0.1. When disk self-gravity is not included, we see again that the disk mass 
has only a minor effect on the disk KL cycle (left column). The oscillation timescales and 
amplitudes are similar. As we discussed in Figured] the reason why the disk mass plays 
a role in the absence of self-gravity is that the disk mass slightly affects the orbit of the 
companion, which in turn affects the onset of the disk KL cycle. With self-gravity included, 
we see a greater suppression of the disk KL cycle by self-gravity than that in Figured! With 
Md = 0.03 M, the KL mechanism barely operates. For the disk mass of = 0.01 M, there 
is a significant reduction. The lowest mass disk remains relatively unaffected by self-gravity. 
Note that for the highest disk mass, the initial disk has a minimum Toomre Q ~ 1.6, 
which puts the disk on the verge of becoming gravitationally unstable to non-axisymmetric 
disturbances. In this paper, we have restricted our attention to simulations in which we do 
not see any sign of gravitational instability, such as fragmentation or clump formation. 


3. Discussion and Conclusions 


Self-gravity introduces a source of disk apsidal precession, in addition to that due to 
the binary companion that causes the KL oscillations. If the self-gravity induced apsidal 


disk precession rate is faster t 
cycle can be suppressed (e.g.. 


ran the binary indu ced precession rate, then the disk KL 


Holman et al. 


19971) . We estimate the magnitude of the 


local self-gravity induced disk apsidal precession rate |5'sg(’")| — 7rGS(r)/(r2(r)r) from 








13 


No Self-gravity With Self-gravity 



Fig. 3.— Similar to Figured] except for different disk aspect ratio H/R = 0.06 and viscosity 


parameter a = 0.1. 
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Equa 


4on (1) of 


Tremaine 


fjl998l) . which is consistent with the rate given in 


Batygin et ah 


(1201 ll) . We determine the binary companion induced disk local precession rate ghi^) 
by approximating the companion’s potential as arising from a uniform ring of mass Mp 
and radius at,, where at is binary separation. Since the disk maintains its coherence, its 
global precession rate is crudely estimated as the angular-momentum-weighted average 
of the local rates (see Equation (9) of Paper II). For the equal mass binary system we 
consider, we take = 0.025a6, rout = 0.3ab (due to disk expansion from the initial 
rout)- The global disk precession rate caused by a binary companion is ~ O.OSflb, 
where fib is the binary orbital frequency. This precession rate is similar in value to the 
nodal precession rate implied by Figure [2] that shows that the disks undergo half the 
nodal precession period at a time of about 15Pb- The global disk precession rate due to 
self-gravity is |^gg| ~ 0.008(Md/0.001M)flb- Requiring |^gg| > |^bl provides an estimate of 
the minimum disk mass to suppress KL oscillations as ~ 0.004M. According to Figure 


7 of 


Batygin et al 


fl201l!) . the critical disk mass for suppressing KL oscillations is generally 


somewhat larger than this value, depending on the disk inclination. Based on that work, 
for the initial disk inclination we have used, the estimated disk mass for KL suppression is 
about a factor of 5 times larger, that is ~ 0.02M. 


In the case of H/R = 0.1, the blue curve in the right column of Figure [U clearly 
shows that the disk KL cycle is weakened for the case of an initial disk mass of 0.03 M. 
In spite of the reduced oscillation amplitude, the disk still undergoes KL cycles. There is 
an eccentricity peak and inclination valley at a time f ~ 7Pb when the disk has a mass 
of about Md ~ 0.02M which we adopt as a crude estimate for the minimum mass for KL 
suppression. 


In the case of H/R = 0.06, the blue curve in the right column of Figure [3] shows that 
the KL cycle is suppressed, for the case of an initial disk mass of 0.03M. In the case of the 
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lower initial disk mass of O.OIM, the KL oscillation appears to be present, but weakened. 
A reasonable estimate of the critical disk mass suppression of KL oscillations is similar to 
our estimate for the other disk model, ~ 0.02M . Therefore, the simulations show that the 
critical disk mass for suppressing KL oscillations is in the range of ~ 0.02M ( 4% of the 
mass of the central star) that is about a factor of 5 larger than the value based on equating 
the self-gravity induced precession r ate with the binary i nduced precession rate. It is in 


agreement with the value implied by 


Batygin et al 


(1201 ll) that is discussed above. 


Papers I and II showed that global disk KL oscillations damp. After damping, the 
disk inclination is at or below critical KL angle. If a planet forms after this stage, then its 
initial orbital tilt is not large enough to trigger the KL mechanism. This reduced tilt poses 
a challenge to the KL model of planets which assumes initially highly inclined planet orbits. 
The tilt evolution of such planet-disk conhgurations requires further study. 


Suppression of the KL cycle by means of self-gravity requires substantial disk masses. 
At the critical disk mass we found, the disk is quite close to being gravitationally unstable. 
In other words, we hnd a relatively narrow window of disk masses where the disk KL 
oscillation is significantly subdued, while the disk is still stable to gravitational instability. 


The highest initial disk mass we consider in this paper, = 0.03 M, shows more 
warping with self-gravity included than in cases without self-gravity. Initially, the 
inclination of the inner parts of the disk increases, while that of the outer parts decreases. 
In Paper II, we found that the critical (minimum) tilt angle for KL oscillations to operate 
in a disk is about 45°, somewhat higher than the critical angle for free particles (39°). In 
this paper, we have limited our analysis to an initial disk tilt of 50° that is just 5° above the 
minimum value for KL oscillations. Such initial conditions produce much more moderate 
KL effects on a disk than occurs at higher tilt angles. We also reported in Paper II that 
strong density enhancements are found during KL oscillations of a non-self-gravitating 
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disk with a larger initial tilt. These density enhancements appear to involve shocks. With 
larger initial tilts than we have considered here, a self-gravitating disk that undergoes 
KL oscillations may undergo disk fragmentation. Fragmentation of a disk is a possible 


mech anism for planet formation and may be aided by binary perturbations (iBoss 


19971. 


20061 ). The disk KL cycle may then provide an alternative means by which a binary 
companion can promote disk fragmentation/clumping. Such effects should be explored in 
the future. 

W.F. and S.H.L. acknowledge support from NASA grant NNX11AK61G. Computing 
resources supporting this work were provided by the institutional computing program at 
Los Alamos National Laboratory. We thank Da niel Price for providing the PHANTOM code for 


SPH simulations and SPLASH code flPrice 


20071 ) for data analysis and rendering of figures. 
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